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ABSTRACT 

Excitation and propagation of waves in a thermally stratified disk with an arbitrary vertical tem¬ 
perature profile are studied. Previous analytical studies of three-dimensional waves had been focused 
on either isothermal or polytropic vertical disk structures. However, at the location in a protoplan¬ 
etary disk where the dominant heating source is stellar irradiation, the temperature gradient may 
become positive in the vertical direction. We extend the analysis to study the effects of the vertical 
temperature structure on the waves that are excited at the Lindblad resonances. For a hotter disk 
atmosphere, the /-mode contributes less to the torque and remains confined near the midplane as it 
propagates away from the resonances. On the other hand, the excitation of the g-modes is stronger. 
As they propagate, they channel to the top of disk atmosphere and their group velocities decrease. The 
differences compared to previous studies may have implications in understanding the wave dynamics 
in a realistic disk structure. 


1. INTRODUCTION 

When an embedded planet orbits around its central 
star, the tidal interactions between the protoplanetary 
disk and the planet will lead to excitation of wav es from 
resonance locations (|Goldreich fc Tremainell97^ . These 
waves propagate radially and may eventually damp or 
dissipate through shocks and transfer angular momen¬ 
tum to the disk. Such pr ocess may also lead to ra¬ 
dial migration of the planet (iLin et alJll996l : iWardllTO^ 
IRafikovI [200211 a nd change of disk structure such as ga p 
formation fe.g.. iTakeuchi et "nillQQOt ICrida et al.ll2006ll . 
Therefore, the properties of these hydrodynamic waves 
and their relation to the disk structure are crucial to the 
understanding of the planet/disk dynamics. 

Since the disks are cold (i.e., the sound speed is 
much less than the circular speed), they are geometri¬ 
cally thin (i.e., the scale height is less than radial dis¬ 
tance to the star, H < r). Yet, the vertical ther¬ 
mal structure is still important to explain the ob- 
served spectral energy d i stribution of the disk s (e.g., 
iChiang fc Goldrei^ll997t iD’Alessio et al.lll998ll . Mea¬ 
surement of such vertical structure in a protoplanetary 
disk has been made possible previou sly with single-dish 
telescopes le. g.. lAkivama et ^1201111 and more recently 
with ALMA ([Rosenfeld et al.ll2013ll . On the other hand, 
recent near-infrared polarization intensity maps suggest 
that the surfaces of some protoplanetary disks are not 
smoot h on the large scale but exhibit concave spiral fea¬ 
tures ([Takami et al.lIMT^ . However, the origins of these 
surface features are still not well-understood as the sur¬ 
face density enhancement requ ired would be too large 
in the 2D linear theories fe.g.. Ijuhasz et all 1201511 . As 
a result, three-dimensional wave theories may shed light 
on the relation be tween the planet-disk interaction and 
observations ('e.g.. lDong et al.lF2015HZhu et al.ll2015ll . 

Three-dimensional wave propagation in accretion disks 
has been studied by several authors in the past. For an¬ 
alytical convenience, these studies assumed a disk that is 
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stable against convection, with either an isothermal verti- 
cal structure (iLubow fc Pringle 1199^ or a polvtropic on e 
([Korvcanskv &: Pringlel Il995t Lubow fc Qgilviel Il998ll . 
Since a polytropic disk contains a finite boundary where 
gas density and pressure vanish above it, such model 
is considered as a wave-guide model due to the bound¬ 
ary conditions applied at the midplane (usually a parity 
condition) and the surface (e.g., vanishing Lagrangian 
pressure perturbation). A set of discrete modes can be 
identified by the vertical eigen function s and the corre¬ 
sponding dispersion relation. iQgilvid (|I998f) provided 
a generalized analysis with magnetic field and a con¬ 
cise classification of waves in a disk, namely, /-mode 
(fundamental), p-modes (acoustic), r-modes (inertial), 
and p-modes (internal gravity). As the waves are ex¬ 
cited at resonance locations wher e the wavelength is 
large ([Goldreich &: Tremaind Il97ll . only /-, g-, and r- 
modes are launched at the Lindblad resonances (LRs) 
where the Doppler- shifted forcing frequency equals the 
epicyclic frequency. ILubow &: Qgilvid (|I998ll studied the 
wave excitation of these modes at the LRs and found 
that most of the angular momentum are indeed carried 
by the /-mode. The authors also found that the wave 
is evanescent in the vertical direction and experiences 
wave-channeling, which is contrary to the previous re¬ 
sults that wa ve refraction oc cu rs due to the grad i ent o f 
sound speed (|Lin et al.l 1199011 . iQgilvie &: Lubowl (|1999f) 
considered a polytropic disk with an isothermal atmo¬ 
sphere. They found that the wave angular momentum 
concentrates near the transition between t he two lay¬ 
ers wh ile propagating radially. Subsequently, iBate et ahl 
(j2002tl studied the nonlinear propagation of axisymmet- 
ric waves using hydrodynamic simulation and verified the 
wave-channeling of these waves. 

While these studies provide a general theory of wave 
excitation and propagation, the cases for a hotter disk 
atmosphere are not explored. In particular, in previous 
models, the temperature decreases along the normal di¬ 
rection to the disk as the disk becomes less optically thick 
and the heat can be easily carried away by radiation. 
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However, at the outer part of the disk, the stellar irradi¬ 
ation may be very efficient to generate a hotter layer as 
the disk surface is flared. 

To study the effects of the vertical temperature gra¬ 
dient, we adopt the aforementioned wave-guide model. 
The background disk structure is assumed to be slowly 
varying in the radial direction (except at LRs) with re¬ 
spect to the waves. Appropriate boundary conditions 
at the midplane and at the top of the disk atmosphere 
are applied to construct a boundary eigenvalue problem, 
with the Doppler-shifted frequency of the wave being the 
eigenvalue. The vertical temperature profile is parame¬ 
terized and modeled such that the temperature varies 
with height from the midplane and gradually reaches a 
steady (isothermal) atmosphere which is heated by the 
stellar irradiation from the central star. The properties 
of the waves, such as wave-channeling, can be studied by 
varying the parameters of the temperature profile. In our 
model, a vertical hydrostatic equilibrium is assumed to 
be consistent with the prescribed temperature profile. As 
we focus on the wave propagation, the change of thermal 
properties of the disk (e.g., variation in the adiabatic in¬ 
dex 7 ) at the top of the atmosphere and the consequence 
of wave dissipation are neglected in this work. 

The paper is organized as follows. In Section [21 we first 
review the formulation and derive the governing equa¬ 
tions and the relevant boundary conditions. In Section 
|3l we present and discuss the result for the wave excita¬ 
tion and propagation. Finally, we summarize this work 
and draw a conclusion in Section |4l 


2.2. Basic State 

The basic state of the disk can be solved by consid¬ 
ering the radial and vertical directions separately. The 
equilibrium velocity is given by Uq = rf 2 (r)e,^, where the 
angular frequency of the gas H(r) is determined by the 
radial force balance at the midplane. 


rH(r)^ 


dr 


-f 

2=0 


1 dPp 
Po dr 


5 

2=0 


( 4 ) 


where po and Pq are the equilibrium density and pressure, 
respectively. The gravitational potential of the central 
star is given by 


$s(r, z) 


GM 

\/r2 -I- ’ 


( 5 ) 


where M is the mass of the central star. For a cold disk, 
the small correction due to the pressure force and the 
gravity of the planet are ignored, so that the equilib¬ 
rium disk is Keplerian (i.e., H oc The vertical 

disk structure is determined by the vertical hydrostatic 
equilibrium. 


dPo 
dz ~ 


( 6 ) 


where = fl{r)'^z is the z-component of the acceleration 
due to stellar gravity in the thin-disk limit. 

Assuming the ideal gas law applies, we have 


2 . FORMULATION OF THE PROBLEM 
2.1. Governing Equations 

A geometrically thin and cold gas disk orbiting around 
a central young stellar object is considered. An embed¬ 
ded planet is treated as a perturber and tidally interacts 
with the disk. Eor simplicity, the self-gravity of the gas 
is neglected although it may be important in the very 
early stage of the disk. The basic state of the system is 
described by an axisymmetric disk in a time-steady equi¬ 
librium without accretion. The inviscid hydrodynamic 
equations reads 

|^+V-(pu) = 0, (1) 

d 1 

"+u-Vu = --VP-V$tot, (2) 

dt p 

where p, u, and P are the gas density, velocity, and pres¬ 
sure, respectively, and $tot is the total gravitational po¬ 
tential. In this work, we use the cylindrical coordinates 
(r, (fi, z) centered at the star. The total potential $tot is a 
sum of the contributions from the central star (subscript 
s) and the orbiting planet (subscript p): 

4>tot = (3) 

where <I>p includes an indirect term due to the choice of 
the origin being away from the center of mass. To close 
the equations, we model the basic state of the disk with 
an empirical formula for the temperature and assume 
the linear perturbation is adiabatic. The effects of the 
magnetic field and viscosity are also ignored in this work. 


Po(r,z) = - po{r,z)To{r,z), (7) 

ma 

where ks is the Boltzmann’s constant, is the average 
molecular mass (in grams), po{r, z) and To{r,z) are the 
equilibrium gas density and temperature, respectively. 
The local scale height H at the midplane is defined in 
terms of the midplane temperature; that is 


= 




ks Tmid 


( 8 ) 


where c and Tmid = T(r, z = 0 ) are the inidplane sound 
speed and temperature, respectively. Eor a cold disk, the 
vertical extent of interest is of order of the scale height, 
but much smaller than the radius. 


z 

r 


H{r) 


e < 1 . 


(9) 


Eor H/r = 0.05 at LRs, the error of gz due to the thin- 
disk approx imat ion is about 0.24 at the top boundary 
(see Section [2^ . 


2.3. Model for Vertical Structure 

Consistent with the assumption of a cold and thin disk, 
we further assume the separation of variables applies such 
that the radial variation of the equilibrium variables (e.g., 
Po, Pq) are determined by their midplane values (e.g., 
Pmid, Tmid) which are labeled with subscript “mid”. The 
midplane variables are varying in a radial length scale 
that is larger than H. Therefore, it is useful to define a 
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with constant temperature, which is given by 


ri + (l/2)A[l-cos^(Z/Za)"] 
\l + A 


A Z<Z, 
ifZ> Z,, 
(14) 


where A is the fractional change of the vertical tempera¬ 
ture, Za is the transition height to the isothermal atmo¬ 
sphere at a temperature T = 1 +A, and n is a parameter 
determining how rapid is the transition. The parameter 
n = 2 is chosen such that the equilium can be solved 
analytically (Appendix [B]) and there is a simple stability 
criterion against convection as we will discuss below. The 
transition from the bottom layer to the isothermal atmo¬ 
sphere is smooth such that T'(Za) = 0 , where the prime 
(and hereafter) denotes the Z derivative. Above Z^,, both 
gas pressure and density fall off as a Gaussian with dis¬ 
tance from the midplane (i.e., oc exp—Z^/(l -f A)). 


Fig. 1. — Vertical profiles of the equilibrium state. The top, 
middle, and bottom panels are the dimensionless profiles for gas 
temperature, pressure, and density, respectively. The solid line 
shows the fiducial model with A = 2.0, Za = 2.0 (indicated by the 
vertical dashed-dotted line). The dashed line shows the polytropic 
disk with s = 5. We denote H as the scale height at the midplane. 


dimensionless 2 coordinate Z = z/H{r), and 


^ Mr,Z) 
TnUr) ’ 

p Pojr.Z) 
Pmid(r) ’ 

^ Pojr, Z) 

^ Pmid(r) ’ 


( 10 ) 

( 11 ) 

( 12 ) 


where T, P, and p are dimensionless counterparts of Tq, 
Pq, and pq, respectively, and are functions of Z. The 
vertical hydrostatic equilibrium in Equation ([5]) can be 
written in the following dimensionless form. 


dP _ PZ 

Iz ~ F’ 


(13) 


where P = pT. Given a vertical temperature profile 
T{Z), the equilibrium pressure and density can be ob¬ 
tained by integrating Equation (1131) and using the ideal 
gas law. 


2.3.1. Temperature Profiles 

In the outer region of the disk, the effect on the temper¬ 
ature due to stellar irradiation dominates viscous heat¬ 
ing. In addition, long-wavelength emissions from embed¬ 
ded large-sized grains and some molecules become op¬ 
tically thin to the disk. Therefore, it is expected that 
the temperature of the disk increases with height from 
the midplane. For simplicity, we do not explicitly solve 
the energy equation with sources of heating and cool¬ 
ing. This allows us to focus on the linear response to 
the disk. Instead, we parameterize the vertical struc¬ 
ture using the relative temperature increase in the disk 
atmosphere, i.e., A = - T'miri)/T'mid- I nspire d by 

iDullemond et al.l l)2002ll and iRosenfeld et all (|20I3ll . our 
temperature model gradually extends into an atmosphere 


2.4. Deviation from Equilibrium 

We proceed to derive a set of equations for the variables 
deviated from the basic state. The velocity perturbation 
can be expressed as Ui = ue^-l-ve^^-l-we^, where u = Uri-, 
V = ityji, and w = Uzi- Keeping the nonlinear terms and 
canceling the terms for the equilibrium disk, the compo¬ 
nent form of Equations can be rewritten as 


dp 

m 


Id Id 

[rpu) -f - 7 ^(pv) 
r or r oip 


dz 


(pw) = 0, 


(15) 


Du 

iDt 

Dv 

~Di 


-2IIv = - 


U'v „ „ 

-h 2Bu = — 

r 


Dw 

~Dt 


I dPi _ d 
p dr dr 
I dPi _ 1_^^ 
pr dip r dip 
i dPi _ 
p dz dz 


where 


_D__a _5 

Dt ~ dt^"^ dr ^ r dip^'^dz 


(16) 

(17) 

(18) 

(19) 


is the Lagrangian derivative, Uip = rfit{r) -|- v is the total 
circular velocity, Pi is the pressure perturbation, 4>p is 
the potential perturbation, and the Oort’s constant B is 
defined by 


r dfl 
^ 2 *^ 


O. 


( 20 ) 


The coefficient 2B in Equation dm) is also commonly 
expressed as K^/2fit, where k is the epicyclic frequency 
defined by = 4BD,. In deriving Equations (fT5l) - (fT8)) . 
no assumption has been made for the amplitude of the 
perturbation. 


2.5. Linear Perturbation 

We proceed to linearize Equations (fT^ - (fT8l) by consid¬ 
ering the linear waves. Assuming the embedded planet 
is orbiting circularly at an angular frequency Lip, the 
planet’s potential can be expressed as a Fourier series 
that is periodic in ip and t, 

^p = ^^mir,z)exp[im{ip-Llpt)], ( 21 ) 
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where m is a positive integer. As Equations (flSjl - lfTSl) do 
not depend on t and ip explicitly, they are linearized and 
Fourier-transformed into the following form 

» 15 / X imonv d , X „ /„„x 

- iujpi +- — {rpQu) ^ -|--^(pow)=0, (22) 

r or r oz 

1 ^ d 

-lWM-2flv-|-- 5 --P 1 = (23) 

po or or 

-ilJjY + 2Bu+—Pi=-—^m, (24) 

por r 

Pi 1 d d 

- lUJW + g—+-—Pi =- — ^rn, (25) 

Po p oz Oz 

where w(r) = m[flp — fl(r)], and the subscripts of the 
Fourier components for pi, u, v, w and Pi are suppressed 
for clarity. The perturbational equation for the specific 
entropy reads 


— iCj 






= 0 , 


(26) 


where 7 is the adiabatic index. 

We focus on the horizontal resonances and ignore the 
vertical forcing due to the gravity of the planet. In the 
thin-disk limit, we ignore the terms due to the azimuthal 
gradient of v and Pi in Equations (1221) and (l24)) , respec¬ 
tively. By eliminating azimuthal velocity v, Equations 
((23l) and (l24)) can be combined into the following equa¬ 
tion, 


(w^ — K^)tt -I- 


ioj dPi 
Po dr 


d^r 


2imQ _ 


dr 




(27) 


2.7. Free Waves 

We first consider the homogeneous solutions (i.e., free 
waves) of linearized equations. These eigenfunctions can 
be then used as the basis functions to expand t he per¬ 
turbational potential (iLubow fc Ogilvielll998ll . Each 
eigenfunction can be categorized into different modes 
by studying the corresponding dispersion relation oj(k) 
where k is the radial wavenumber. As the torque is pro¬ 
portional to the fraction of each mode can be com¬ 
puted. To take into the account the linear term of x 
in Equation (1281) in the proximity of the LRs, we make 
use of the following sepa ration of variables suggested by 
ILubow fc Ogilvid (|i998[ ). 


u{x, 

z) 

= Ni{qx)u{z), 

(30) 

v(x. 

z) 

= iA\(qx)'v{z), 

(31) 

w(x. 

z) 

= qN\{qx)w[z), 

(32) 

Pl{x, 

z) 

= iqk\ {qx)Pi{z), 

(33) 


z) 

= iqk\ {qx)pi{z), 

(34) 


where x is the dimensionless distance from the LR, q is 
the effective radial wavenumber, and Ai(a;) is the Airy 
function satisfying the following ODE, 


Ai"(x) — xAi(a;) = 0. (35) 

In the absence of forcing, Equations (l22]) - (l26)l can be 
combined and simplified into two ODEs by using Equa¬ 
tions O-O- We define the following dimensionless 
variables 


where the forcing on the right-hand-side of the equation 
is only dominant at the LRs. 


2.6. Local Approximation and Normal Modes 

In order to study the vertical wave structure, we con¬ 
sider a radially restricted region such that the radial 
variation of the equilibrium state (subscript 0 ) can be 
neglected. This local appro ximation, or equivalently 
shearing-sheet approx imation (|Goldreich &: Lvnden-Bell 
119651 l^odmanl1199311 . allows the separation of variables 
in the radial and vertical directions in the governing 
equations. In particular, we assume the equilibrium state 
is a function of z only and obtain a set of linear ordi¬ 
nary differential equations (ODEs) in z. We consider the 
waves excited at the LRs, where w(r = vl) = ±k. 


2.6.1. Near the Lindblad Resonances 

In the proximity of the LRs, the wavelength of the per¬ 
turbation is comparable to radial extent of the forcing. 
Therefore, the usua l plane wave approxima tion breaks 
down and we follow iLubow fc Ogilvid ()1998ll to expand 
the singular term to the linear order of dimensionless 
displacement a; = (r — r’i,)/ri, that is 

~ Six, (28) 


where S = rd/dr{K^ — w^)|ri. For a Keplerian disk, 
S = —3(1 ± m)Lt‘j, where {1l = 0(r = r^) and m > 0. 
Thus, Equation (1^ becomes 


„ iuj dPi „ 

— Sxu -\ -— = —lUJ 


2imLl . 


Po dr 


dr 




(29) 


X{Z) = ujPi{z)/{n^HHo), 

(36) 

W{Z) = w(z)/(DiJ), 

(37) 

F = uj/VL, 

(38) 

Q = {Hlrfl\, 

(39) 

S = SIVl^. 

(40) 


After some algebra and making use of dimensionless vari¬ 
ables, we obtain the following differential equations, 


dX 

Hz 

dW 

Hz 


Ar2 


x+[f‘^ - iV^] W, 


A-^ 

7 P] 


P' 

X-—W, 

jP 


(41) 

(42) 


where = 1 at LRs, A = Q^/S is a unknown parameter, 
and 


iV 2 


P_ 

"f dZ P '1’ 


(43) 


is the square of the dimensionless Brunt-Vaisala fre¬ 
quency (normalized by fl). The equilibrium variables 
are obtained by Equation (TT^ in Section 12.31 When 
the problem is solved with appropriate boundary con¬ 
ditions, the parameter A is determined as an eigen- 
value. Equations (HTl) an d (1421) follow those obtained by 
ILubow &: Ogilvid (1199811 except that we assume an ar¬ 
bitrary hydrostatic background profile instead of a poly¬ 
tropic disk. We rewrite the equations in a form that does 
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not require solving Equation (US in advance, i.e., 


dX 

Iz 


+ 

dW _ 

Hz ~ 


- 


7-1 


7 

7 — 1 \ Z 


7 


A- 


1 


T 

X 


T 

T' 

Y 

Z 


ZT' 

~T~ 


W 


X, 


iT 




(44) 

(45) 


where we have used 





T 



(46) 


For non-decreasing temperature profiles {A > 0) in Equa¬ 
tion (fT4|l {T'{Z) > 0), the disk is stable against con¬ 
vection > 0. For A < 0 and n = 2 in Equa¬ 
tion (HI, the stability criterion against convection is 
|4l| < (7 — l)Za/7'^ 5 for our choice of parameters 

n = 2 and Za = 2. Thus, the background is stable 
against convection in our explored parameter space. 


2.8. Boundary Conditions 

Equations (HH) - (H5]) are the governing differential equa¬ 
tions near the LRs. Away from the LRs {F ^ 1), the 
equations can be obtained by replacing A with (F^ — 
K^) (see Appendix lAl Equations (lAll) and (EU). They 
are numerically integrated and are subjected to boundary 
conditions at the midplane (Z = 0 ) and at an arbitrary 
defined surface (at Z = Zb). As the system is symmet¬ 
ric in Z in our simple model, the waves take one of the 
parities. We consider only the even solution where 


dX 

Hz 


z=o 


0 and W{Z = 0) = 0, 


(47) 


in which the pressure perturbation is an even function in 
Z and the warping modes are ignored. 

The adiabatic index is assumed to be constant every¬ 
where (7 = 5/3), such that both the sound speed and 
buoyancy frequency are continuous across different lay¬ 
ers. The upper boundary condition is applied such that 
the Lagrangian perturbation of pressure goes to zero at 
the top of the atmosphere Zb. In terms of the dimen¬ 
sionless variables, it reads 


Ab - WbZb = 0, (48) 


to study the behavior of each mode, we apply a upper 
boundary condition well above the transition height, such 
that the /-mode and g-modes appear to be discrete. For 
r-modes, its spe ctrum remains discrete a s in the locally 
isothermal disk (ILubow &: Pringf^ 1199311 . In this work, 
we follow the convention to label each mode of the same 
type by the number of nodes in the pressure perturbation 
(i.e., X). 

The background equilibrium density and pressure are 
obtained by solving the equations for hydrostatic equilib¬ 
rium. The basic vertical structure is generally smoother 
than the polytropic disk counterpart (see Figured]). A 
reference model for A = 2.0 and Za = 2.0 is used, as 
shown in Figure [2| The parameters are varied in order 
to understand their effects on wave excitation and prop¬ 
agation. 



0.5 1.0 1.5 2.0 2.5 3.0 

r/ro 


Fig. 2.— Temperature profile for A = 2.0. The horizontal and 
vertical axes are r and 2 : in the unit of corotation radius vq. The 
color scale shows the relative temperature with respect to midplane 
temperature Tq at r = tq. The dashed and solid gray lines are the 
transition location {Za = 2.0) and upper computational boundary 
{Zfy = 8.0), respectively. The vertical black dashed lines denote 
the LRs for m = 2. Here we adopt a power-law T^[^{R) oc 
and the aspect ratio H/r = 0.05 at the outer LR. 


At ^ = 0, our model reduces to the locally isother¬ 
mal disk with adiabatic response. This was previously 
studied by ILubow &: Pringld (1199311 , in which a bound¬ 
ary condition is applied such that the wave energy den¬ 
sity is finite at 2 ; —>■ oo. In our case, since an artificial 
upper boundary applied, the discrete g-modes are recov¬ 
ered. However, g-modes are only weakly excited in this 
case. 


where Ab and Wb are evaluated at Z = Zb. 

3. RESULTS AND DISCUSSION 

A local model is constructed by assuming slow radial 
variation of the equilibrium state. For simplicity, the ver¬ 
tical equilibrium temperature profile is taken as a smooth 
continuous function. The functional form is chosen to 
mimic the vertical struct ure determined by the radia¬ 
tive transfer model (e.g., iDullemond et al.l 1200211 . The 
temperature profile is characterized by three parameters, 
namely, A the fractional increase of the temperature of 
the atmosphere relative to that of the midplane, Za the 
transition height to the atmosphere, and Zb the upper 
boundary. In principle, the isothermal atmosphere ex¬ 
tend to the infinity, in which the internal gravity (g) 
modes have a continuous spectrum. However, in order 


3.1. Torque Fraction 

At the LRs, the torque due to the external poten¬ 
tial can be evaluated as (iGoldreich fc Tremaind [l97^ 
ILubow &: Ogilvid 1 199811 . 


= —sgn(w) 


TT^rnEL'I'm 


(49) 


where Sl 


f pdz is the column density at the LR, and 


= r- 


d^n 


dr 


2mfl 


(jj 


(50) 


The formula from two-dimensional analysis is valid here 
because the vertical forcing of the planet is ignored. By 
projecting onto the basis eigenfunctions, the fraction 
of the torque due to each mode of waves can be obtained. 
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W e denote fj the fraction of torque carried by each mode 
j (|Lubow fc Ogilviall998li . which is given by 


fj 


PQ^CjdZ 


/El / poX^dZ, 


(51) 


in the torque integral and results in a reverse trend of fj 
with Zb. For ri-mode, the torque fraction is independent 
of Zb- This is because X is flat and close to zero for large 
Z and so is not sensitive to Zb- 


where we used the relation between u and dimensionless 
pressure perturbation X: 



which is obtained from Equations (l2^ and (1361) . 

3.1.1. Dependence on Temperature Gradient 

For a simpler analysis, the dimensionless transition 
height Za is kept fixed at 2.0. Therefore, the temper¬ 
ature gradient in the layer near the midplane is propor¬ 
tional to the temperature parameter A. The parameter 
changes both the background density po and eigenfunc¬ 
tions {X,W). As indicated in Equation (ISTI) . the torque 
fraction is determined by the overlap integral between po 
and X. Since fj = 1 among all modes, the torque 
fraction of any particular mode also depends on that of 
others. Most of the torque is carried by the first few 
lower-order modes (i.e., /, ri and pi, where the subscript 
1 means one node in X), although other higher-order 
modes are also excited when A is large. Therefore, we 
focus on these lower-order modes for the analysis. The 
dependence of fj on A for different values of Zb is shown 
in Figure [3l 

For /-mode, the torque fraction increases with decreas¬ 
ing A and becomes less sensitive to Zb- For pi mode, 
the torque fraction peaks at A ^ 0.5 and is generally 
larger for larger Zb- On the other hand, the ri mode 
contributes more torque than the pi-mode when A < 0 . 
While its torque fraction also peaks at some finite value 
of A, the torque fraction is independent of Zb (Figure [3]). 

The general increasing trend of // for decreasing A 
(< 1 ) can be explained by the decrease of fj for pi- and 
ri-modes. In particular, po becomes more concentrated 
near the midplane and so the torque integral in Equation 
m depends sensitively on the overlapping of po and 
X for small Z. Either a smaller numerical value (e.g., 
pi-mode) or an oscillatory solution (e.g., ri-mode) can 
result in a small torque fraction. 

3.1.2. Dependence on Upper Boundary 

Since po is independent of the location of the upper 
boundary Z}, in our model, the torque fraction depends 
only on the functional form of X. For /- and pi-modes, 
the eigenfunctions X are similar to those in a polytropic 
disk, where X{Z) is generally an increasing function. 
While X{Zb) is a local maximum, the value of X for 
/-mode at the midplane varies (without the effect of am¬ 
plitude) and becomes very large at some frequency F. 
This eventually relates to the wave-channeling that we 
will discuss in the next section. Note that X has zero 
and one node for /- and p-modes, respectively. Larger 
Zb generally results in a smaller overlapping of X and po 
near the midplane. Thus, fj decreases with increasing Zb 
for /-mode and pi-mode (for small A). However, when 
A is large, the amplitude of oscillation near the node for 
pi-mode is also large. This cancels out the contribution 


3.2. Wave Propagation 

The wave propagation of each mode is determined 
by the corresponding dimensionless dispersion relation 
F{K), where F = uj/Ul and K = kFd. In Figure SI 
the dispersion relations are shown for /-, pi-, and ri- 
modes {Z^, = 2 and Zb = 8). The slope of F(K) for 
/-mode increases with larger A. It is expected for it s 
acoustic character for AT —>• 0 (iLubow fc OgilviellI998f) . 
although there is no simple WKB relation for the mod¬ 
erate values of K. On the other hand, as p-modes are 
driven by buoyancy (|Ogilvielll998fl , the forcing frequency 
approaches the Brunt-Vaisala frequency near the top 
of the atmosphere for moderate values of K such that 
F{K) ^ N (X 1/(1 -I- A)^/^ (see Equation (l43l) l. Using 
the temperature in the atmosphere, the ratio between F 
for A = —0.5 and 2.0 is about ~ 2.4, which is similar 
to the estimate from Figure SI 

As the Doppler-shifted frequency w is a function of r 
for a given azimuthal wavenumber to, the radius can be 
expressed as 

L = + (53) 

ro 

where rg is corotation radius. Note that F is positive for 
r > ro and negative for r < rg. The LRs are located at 
r^/rc = (1 ± 1 /to)^/^, where the plus and minus signs 
indicate the outer and inner LRs, respectively. For the 
purpose for illustrating wave propagation, we consider 
only the outer LR [tl = r^). In Figure SI we show 
the frequencies for /- and pi-modes that propagate from 
the outer LR with to = 2 (note that r-modes are absent 
in this region). In particular, the wavenumber pi mode 
rapidly increases as it propagates outward. As the group 
velocity of pi-mode decreases significantly, the wave am¬ 
plifies at some finite distance from the LRs (r < 2rb) 
due to the wave action conservation explained below. In 
particular, the location where gi-modes “piled-up” are 
closer to the LRs for the increasing A. For realistic disk 
with positive A and higher to, we expect such effect oc¬ 
curs closer to the LRs. 

For waves that propagates in a non-uniform moving 
medium, the conservation of ener gy or angular momen¬ 
tum wave action can be derived (|Bretherton fc GarrettI 
119681 [Andrews fc Mclntvrelll97j^ . given that the mean 
flow can be meaningfully defined. A similar application 
regardin g spiral density wav es iii galac tic disk was dis¬ 
cussed in lToomr'i (|196^ and iShul ( 197011 . In the absence 
of self-gravity, we only need to consider the flux of wave 
action due to advective transport. In any case, as the ra¬ 
dial flux of angular momentum wave action is conserved 
in the absence of dissipation, it can be used to relate 
the amplitude of the waves at different radii in the lin¬ 
ear theory. For single-mode WKB (radial) waves (e.g., 
Pi oc expf J k(r)dr), the time-averaged, vertically inte¬ 
grated flux of angular momentum wave action can be 
derived from Equations (IMl) and (|24l) and is given by 
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AAA 


Fig. 3.— Torque fractions of the /-mode (left), pi-mode (middle), and ri-mode (right) versus temperature parameter A for different 
values of the height of upper boundary The vertical dashed line indicates the locally isothermal disk (A = 0). The transition height is 
at Za = 2.0. Note that for ri mode, the torque fractions are the same for different Zb- 




Fig. 4.— Dispersion relation of different modes (Za = 2,Zb = 8) 
are shown for A = —0.5 and 2.0. The solid, dashed, and dotted 
lines represent the /-, gi- and ri-modes, respectively. The waves 
are launched at LRs where F = 1 and A' = 0. 


(e.g., iLubow k, OgilvidllQQSD 


p(.a-) — 


irrm f 


pH'^dz, 


(54) 


where k is the radial wavenumber. Note that this equa¬ 
tion is general for WKB theories with k'^ ^ m?/r"^ (i.e., 
valid away from LRs) and is derived without our previ¬ 
ous assumption in Section [2.5l to neglect azimuthal pres¬ 
sure gradient when computing the eigenfunctions. By 
expressing the above formula in terms of dimensionless 
variables, we have 


= TTTO 


^ ir2(^2_^2) 



p{Z)X'^{Z)dZ, 

(55) 


where Pmid(?') is the midplane gas density and k = 1 is 
the normalized epicyclic frequency for a Keplerian disk. 
Here we denote a = a(r) as the slowly-varying ampli¬ 
tude of the waves (i.e., u{r,z) = au{z)e'^’^'~ at non-LR 
locations). On the left-hand-side of Equation (1551) . the 


Fig. 5.— Plot of radial wavenumber versus radius for A = —0.5, 
0.0, and 2.0 (Za = 2,Zb = 8). The solid and dashed lines represent 
the /- and gii-modes, respectively. 


value of for each mod e is determined by its corre- 
spending value at the LRs (iLubow fc Ogilvielll998D . i.e., 

where and fj are the fraction and total 
value of the torque in Equations and (ED, respec¬ 
tively. 

The confinement and wave-channeling of the waves are 
examined by considering the density of angular momen¬ 
tum wave action, which is given by 

^^“^ = ^P(H^ + Iwp), (56) 

where averaging over one wave period is assumed. Conse¬ 
quently, the vertically-integrated flux £/^‘^^dz, 

is related to the radial flux by pA) = A^^^Vg, where 
Vg = dCj/dk is the radial group velocity. For a monochro¬ 
matic wave (w/m = Hp being a constant), the relevant 
conservation law for the wave action reads 


dr 


= 0 , 


(57) 


which implies being a constant except at the LRs. 


















































For a given disk model and a dispersion relation of a 
particular mode, the wave amplitude can be expressed as 
a function of r, which is obtained by keeping a con¬ 
stant and using normalization such as f p^X^dZ = 1 at 
each radius. Here, we use simple power laws for midplane 
density (pmidir) oc r~P) and temperature {Tmidir) oc 
r“'J). In the following, p = 2.25 and q = 1/2 such that 
the surface density goes as Sq oc pqH oc r~^. In Figure[6l 
the density of wave action of /-mode are shown on 
the r-z plane for A = —0.5, 0.5, and 2.0 and in a radial 
range oi < r < 2rL. When the atmosphere is colder 
than the midplane {A < 0), the density of wave action 
peaks at the midplane and gradually channels near the 
transition where Za = 2.0. This wave- c hann eling effect 
was also reported bv iQgilvie fc Lubo^ (|1999ll where an 
isothermal atmosphere is attached to a polytropic disk 
near the midplane. Since the sound speed is continu¬ 
ous in our model across Za, the density of wave action 
do not stay at a particular height as in their work. For 
the cases with a hotter atmosphere {A > 0), the den¬ 
sity of wave action peaks near (but not exactly at) the 
midplane. Instead of channeling to the base of hot at¬ 
mosphere in the disk, the /-mode channels towards the 
midplane where the sound speed is lowest. This sug¬ 
gests that the direction of wave-channeling for /-mode 
depends on the temperature gradient, which is expected 
from its acoustic nature. For gi-mode shown in Figure 
[7l channels even higher as it propagates away from 
the LR and peaks close to the upper boundary. 

While the vertical velocity perturbation w{z) peaks at 
the upper boundary for /- and g-modes, the peak height 
of the density of angular momentum wave action dif¬ 
fers. The concentration of wave energy locates at differ¬ 
ent heights due to the wave-channeling. Although the 
nonlinear effect is not included in the current work, the 
wave action density shown in Figures [6] and [7] indicates 
the differences among different modes. 

The upper boundary serves as a lid which allows sur¬ 
face modes to exist. This is only an idealization such 
that the waveform does not extend to inhnity. This is 
expected that, far above the midplane, the thin-disk ap¬ 
proximation breaks down and the thermal property of the 
disk changes. At a greater height, where the heat from 
stellar irradiation can be easily radiated away, the ther¬ 
mal time is expected to be so short that the gas response 
be comes isothermal. This case was previously studied 
bv IQgilvie fc Lubo^ (|I999ll where the isothermal gas re¬ 
sponse (7 = 1 ) in the atmosphere is considered and a 
matching condition can be derived at the base of this 
layer. When such boundary condition is adopted (with a 
corresponding matching condition), no g-mode is found 
using our method. This is consistent with the specula¬ 
tion that 5 -modes are in the continuous spectrum. 

Among the lower-order modes considered in this work, 
the eigenfunctions are mostly evanescent in the vertical 
direction. The vertical wavenumber is thus not well- 
defined. The vertical group velocity is essentially zero 
even though the number of nodes may change slightly 
(e.g., from zero to one) during propagation. The pres¬ 
ence of wave-channeling instead of wave propagation 
in the vertical direction (i.e., refraction) was also con- 
firmed in the asym metrical hydrodynamic simulations in 
iBate et ^ (|2002f ) for a polytropic disk. 


4. SUMMARY AND CONCLUSION 

In this work we explore the properties of the waves 
that are excited at the LRs in a planet-disk system. In 
particular, the disk models with a hotter atmosphere is 
studied. The linear theory under the WKB approxima¬ 
tion is developed based on the earlier work by Lubow and 
Pringle and the separation of variables in terms of radial 
and vertical coordinates is applied. This model is essen¬ 
tially an 1-l-lD calculation, in which the eigenfunctions 
are calculated at each radius r (or equivalently at each 
Cj) and are connected depending on the radial structure 
of the disk by means of the conservation of wave action 
(see. Equation [57]) . 

At each radius (including LRs), the local boundary 
value problem is solved for the eigenvalues and eigenfunc¬ 
tions. Each eigenvalue corresponds to different mode, 
which can be categorized ba sed on their behavior for 
large radial wavenumber (k) (|Ogilviel|l998ll . For exam¬ 
ple, the linear dispersion relation of /-mode for large 
k reflects its acoustic nature. In Section 13.11 the eigen¬ 
functions are used to compute the torque fraction of each 
mode. Unless A is very large, the /-mode carries most 
of the angular momentum flux away. The temperature 
of the disk atmosphere and its gradient affect the back¬ 
ground density and sound speed, which in turns change 
the contribution of different waves excited at LRs. As the 
wave propagates away from the LRs, the density of the 
wave action concentrates at different height. This wave¬ 
channeling phenomenon in the vertical direction differs 
from the usual wave propagation as the group velocity is 
only tangential to the disk plane (i.e., dCo/dkz = 0 ). 

For a disk atmosphere hotter than the midplane {A > 
0 ), we hnd that the /-mode remains confined within the 
layer closer to the midplane (Figure El), whereas the g- 
modes retain its surface-wave character and continue to 
channel upward (Figure [7]). Since the group velocity of 
/-mode is almost constant when it propagates far enough 
from the LR, the concentration of is mainly due to 
the change in eigenfunctions. On the other hand, as a 
result of the angular momentum conservation and the 
decrease in group velocity (Figured]), the amplitudes of 
the g- and r-modes increase as they propagate. Depend¬ 
ing on the strength of the tidal forcing (i.e., mass of the 
planet), the increasing amplitude may lead to nonlin¬ 
ear dissipation. Combining the effect of wave-channeling, 
our understanding on the formation and structure of the 
spiral shock may change for a three-dimensional disk. 

We parameterize the vertical temperature prohle T(Z) 
for an equilibrium disk using Equation (1141) . in order to 
mimic the re sults from previous rad iative transfer calcu¬ 
lations le.g.. iDullemond et al.ll200^ . By expressing the 
governing equations in terms of the arbitrary tempera¬ 
ture profile T(Z), one can easily apply our formulation to 
study different models used in observations and simula¬ 
tions. However, identification of different modes for each 
m in a realistic planet-disk system remains a challenge 
and will be left for future investigation. In conclusion, 
we study the behaviors of different waves that are excited 
at LRs. The properties of these linear waves will help us 
better understand the tidal interactions between planet 
and disk. 


4.1. Application to Nonlinear Effects 
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r/ri r/r£ r/ri 


Fig. 6 .— Variations of density of angular momentum wave action for /-mode. The linear color scale in the arbitrary on each plot. The 
values of A are —0.5, 0.5, and 2.0, for the left, middle, and right figures, respectively. The horizontal axis is between rj^ and 2r£^. 



rjrL r/ri 


Fig. 7.— Same as Figure [ 6 ] but for 5 ' 1 -mode. 



Lastly, the linear study of waves may still give us some 
hints about the nonlinear processes. There are a few 
factors determining whether the waves will eventually 
break, for example, the strength of the forcing, local 
sound speed, and the character of the mode. For /- 
mode, the wave energy are confined within the central 
colder layer, but the amplifications are only moderate 
and not very sensitive to the choice of A (see Figure H]). 
For 5 -modes, which carry more angular momentum in 
the case of hotter disk atmosphere, they amplifies signif¬ 
icantly and channel into the top of disk atmosphere. It 
will be interesting to investigate the nonlinear behavior of 
5 -modes in the atmospher e. Other dissipation p rocesses 
such as viscous d amping ([Takeuchi et al.|j l99fiH and ra¬ 
diative damping (jCassen fc Wooluml 119961 1 may be rele¬ 
vant to the linear waves in this work as well. 

Together with the finding that the torque fraction car¬ 
ried by each mode depends sensitively on the tempera¬ 
ture gradient, we expect that the shock structure (e.g., 


shock location) in such a three-dimensional disk may 
be very different from the two-dimensional disk. Con¬ 
sequently, the change in the location of the deposit of 
angular momentum may lead to differences in terms of 
planet migration as well. In particular, the excitation 
of 5 -modes in a thermal-stratified disk may modify what 
we expect from an isothermal disk in the literature (e.g., 
iTanal^ et al.ll200^ . 

Nevertheless, recent three-dimensio nal numerical sim - 
ulations of planet-disk interaction by iZhn et rdl ()2015r i 
show that the inner spiral shock (within corotation ra¬ 
dius) are formed at a few scale height from the midplane 
(but not the computational boundary). While the wave¬ 
channeling effect in an initially-isothermal disk adopted 
in their study is expected to be insignificant, the non¬ 
linear effects from the shock-heating at the surface is 
unknown. Therefore, our study on the wave channeling 
may provide a basis for future investigations. 

The work is supported by a MOST grant in Taiwan 
through MOST 103-2I12-M-001-027. 


APPENDIX 

A. AWAY FROM THE LINDBLAD RESONANCES 

Away from the (horizontal) resonance locations where the waves are excited, the radial wavelength of the perturbation 
is small compared to the variation length scale of the background equilibrium. The pert urbational equations are very 
similar to the case nea r LRs, except there is no explicit x dependence in Equation (l?7l) . iLubow &: Pringle! (II993I1 and 
iKorvcanskv fc Pringle! (j 19951 1 performed the local calculations of the axisymmetric (m = 0) free waves in the vertically 
isothermal and polytropic disks, respectively. In order to perform the normal mode analysis for non-axisym metric 
wayes, we further neglect the adyection terms due to the background shear as done in ILubow fc Ogil\^ (|1998[ 1. One 
reason is that the mixture of t ime (t) and space (r) does not allow the normal modes to exist in the form of shearing 
wayes (|Ryu fc Goodm^ll992ll . Since we are more interested in the yertical structure of wayes at different radii, the 
information on the radial front is sacrificed. In this limit of WKB approximation in the radial direction, we apply the 
Fourier transformation in r for the perturbational yariables (i.e., expifcr). By keeping (li{r) as a local constant, we 
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Fig. 8 .— Variations of density of angular momentum wave action for /-mode, in physical space. The linear color scale in the 

arbitrary on each plot. The values of A are —0.5 and 2.0, for the left and right panels, respectively. The blue dotted line indicates the 
transition height at z = 2H(r). The horizontal axis is between and 2rL for m = 2. The corotation radius is at 100 AU. The power-law 
index for density and temperature are p = 2.25 and q = 0.5, respectively. 


effectively dropped the differential shear in the disk (I Vishniac fc DiamondlflOS^ . Therefore, using the dimensionless 
variables, the differential equations become 


where we define 


dX 

Hz 

dW 

Hz 


—X + {F^ 


N^)W, 


■ ^2 
F 2 - 


PjZ) ' 
lP{Z)_ 


X- 


P'jZ) 

lP{Z) 


w, 


X(Z) 


iCjPi{z) 

WhHHh 


W{Z) 


w{z) 
flH ’ 


K = kH, 


(Al) 

(A2) 


(A3) 

(A4) 


and H is the scale height at the inidplane (Equation ([S])). In terms of T{Z), we have 


dX 

[/7-A 

z , r'(z)' 

dZ ~ 

.V 7 j 

T(Z) ' 

r(z)J 

dW 


I 

A + 

dZ 

_E2 - k2 

7T(Z). 


A 


lT{Z) 


- 


W. 


7-l\ Z^ ZT'(Z) 


7 J T(Z) T(Z) 


W, 


(A5) 

(A6) 


B. VERTICAL HYDROSTATIC EQUILIBRIUM 

For n = 2 in Equation (I14p . the vertical hydrostatic equilibrium can be solved analytically. The dimensionless 
equilibrium pressure is given by 


P(Z) = 



: arctan 


2FT+A 

zl z'^-Z 


\/\ + AtanTT 


if Z<Z, 


ifZ> Za, 


' 2vTT:4 2(1+A) 

where Z^ is the transition height. The equilibrium density is given by p(Z) = P(Z)/T(Z). 


(Bl) 
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